Efficiency of Nonlinear Particle Acceleration at Cosmic Structure 

Shocks 



Hyesung Kang 

Department of Earth Sciences, Pusan National University, Pusan 609-735, Korea 

^ ■ kang@uju . es . pusan . ac . kr 

O ' 

<N ; and 
-i— > 
O 

Q ■ T.W. Jones 

^ ■ Department of Astronomy, University of Minnesota, Minneapolis, MN 55455 

twj @msi . umn . edu 

(N 

' ABSTRACT 

O 

O . We have calculated the evolution of cosmic ray (CR) modified astrophysi- 

cal shocks for a wide range of shock Mach numbers and shock speeds through 

6 

well as pre-existing, upstream CR populations. Bohm-like diffusion is assumed. 
We model shocks similar to those expected around cosmic structure pancakes as 



numerical simulations of diffusive shock acceleration (DSA) in ID quasi-parallel 
plane shocks. The simulations include thermal leakage injection of seed CRs, as 



well as other accretion shocks driven by flows with upstream gas temperatures in 
the range T = 10 4 — 10 7,6 K and shock Mach numbers spanning M s = 2.4 — 133. 
We show that CR modified shocks evolve to time-asymptotic states by the time 
injected particles are accelerated to moderately relativistic energies (p/mc ^ 1), 
and that two shocks with the same Mach number, but with different shock speeds, 
evolve qualitatively similarly when the results are presented in terms of a char- 
acteristic diffusion length and diffusion time. We determine and compare the 
"efficiencies" of CR acceleration in our simulated shocks by calculating the ratio 
of the total CR energy generated at the shock to the total kinetic energy that 
would pass through the shock over time in its initial frame of reference. For these 
models the time asymptotic value for this efficiency ratio is controlled mainly 
by shock Mach number, as expected from the aforementioned similarity in CR 
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shocks. In the presence of a pre-existing CR population, shock evolution proceeds 
similarly to that for higher thermal injection rates compared to thermal leakage 
CR sources alone. This added contribution has little or no impact on the post- 
shock or CR properties of the high Mach number shocks simulated. The modeled 
high Mach number shocks all evolve towards efficiencies ~ 50%, regardless of the 
upstream CR pressure. On the other hand, the upstream CR pressure increases 
the overall CR energy in moderate strength shocks (M s ~ a few), since it is a 
significant fraction of the shock ram pressure. These shocks have been shown to 
dominate dissipation during cosmic structure formation, so such enhanced effi- 
ciency could significantly increase their potential importance as sources of cosmic 
rays. 

Subject headings: acceleration of particles-cosmic rays- hydrodynamics- meth- 
ods:numerical 



1. Introduction 

Collisionless shocks form ubiquitously in tenuous cosmic plasmas via collective, electro- 
magnetic viscosities. The formation process of such shocks inevitably produces suprathermal 
particles in addition to thermal particles with Maxwellian velocity distributions (Blandford 
and Eichler 1987; Jones and Ellison 1991; Draine and McKee 1993). These nonthermal parti- 
cles can be further accelerated to very high energies through the interactions with resonantly 
scattering Alfven waves in the converging flows across a shock, i.e., by diffusive shock accel- 
eration (DSA) (Drury 1983; Blandford and Eichler 1987; Malkov and Drury 2001). Detailed 
nonlinear treatments of DSA predict that a small fraction of incoming thermal particles can 
be injected into the CR population, and that a significant fraction of the shock kinetic energy 
can be transferred to CRs (e.g., Ellison, Baring & Jones 1995; Kang, Jones, and Gieseler 
2002). 

Simulations of DSA in spherical supernova remnants (SNRs) indicate that CRs can 
absorb up to 50% of the initial blast energies (e.g., Berezhko and Volk 1997, 2000). Support 
for rapid and efficient DSA in that setting has been provided by recent X-ray observations 
of young SNRs such as SN1006 and Cas A that indicate the presence of short-lived, TeV 
electrons emitting nonthermal synchrotron emission immediately inside the outer SNR shock 
(Koyama et al. 1995; Allen et al. 1997; Bamba et al. 2003). On a galactic scale it is well 
known that the CR energy density is comparable to the gas thermal energy density in the 
interstellar medium and plays important dynamical roles in the evolution of our Galaxy. 
Although the Galactic CRs are commonly believed to be accelerated mostly at SNR shocks, 
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CR acceleration is probably important in all shock-heated astrophysical plasmas. 

CR populations are also indicated by extended, diffuse nonthermal emissions in some 
galaxy clusters (e.g., Miniati et al. 2001, and references therein). One likely contribution 
to that CR population is large scale shocks associated with the formation of the clusters. 
According to hydrodynamic simulations of large scale structure formation in the universe, 
nonlinear structures such as pancakes, filaments, and knots are surrounded by "external" 
accretion shocks and contain complex webs of "internal" flow shocks including, but not 
limited to, accretion shocks around individual clusters and merger shocks inside clusters 
(Miniati et al. 2000). Recently, Ryu et al. (2003) studied the characteristics of such cosmic 
shocks and showed that they have a wide range of physical parameters with shock speeds, 
u S) up to ~ 3000 km s -1 , preshock temperatures of 10 4 < T < 10 8 K, and Mach numbers up 
to a few 100. They showed that in the present universe the mass inside nonlinear structures 
has been shocked approximately twice on average over cosmic time, and that shocks with 
2 ^$ M < 4 have contributed ~ 1/2 of the total energy dissipated at the shocks. Utilizing 
nonlinear DSA model calculations of Kang, Jones, and Gieseler (2002), Ryu et al. estimated 
that the ratio of CR energy to gas thermal energy resulting from dissipation at cosmological 
shocks is ~ 1/2. Hence, the CR energy density could be dynamically important in the 
intracluster medium, just as in the interstellar medium of our Galaxy and might have had 
some dynamical influences on large scale structure formation. 

The above "cosmic structure shocks" are an inevitable consequence of gravitational col- 
lapse, and are central to virialization of diffuse baryonic matter. The DSA process behind 
our present discussion depends, of course, on the existence of a magnetic field in the shock 
vicinity to mediate Alfvenic turbulence that scatters the CRs. The effectiveness of the scat- 
tering at a particular particle momentum, p = •yvm, can be expressed in terms of the DSA 
diffusion time, t d (p) = n(p)/u 2 s , where k = (l/3)r g v(B /5B) 2 is the spatial diffusion coeffi- 
cient, r g = pc/(eB) is the particle gyroradius, B is the magnetic field, and 5B characterizes 
the magnitude of resonant wave magnetic fields. Particles will generally be accelerated from 
low energies to E ~ pc on timescales a few times td(p) (see equation 2). Efficient transfer 
of kinetic energy to CRs typically becomes apparent in DSA simulations by the time tran- 
srelativistic CRs are accelerated, establishing the relevant time scale for the introduction of 
nonlinear feedback into the shock dynamics. 

Evidence is strong that magnetic fields exist in clusters at levels exceeding at least a 
few xlO~ 7 Gauss (e.g., Kronberg 2001). Fields of that strength lead, with Bohm diffusion 
(SB = B), to diffusion times for transrelativistic CRs, t d = (l/3)pcv/ (eBu 2 ) of the order of 
a few years for p ~ mc in fast cosmic shocks with speeds ~ 10 3 km s _1 . 

Direct evidence for microGauss-level magnetic fields where most cosmological shocks 
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form, that is, outside individual clusters, is currently quite limited. However, there are good 
reasons to expect fields adequate for efficient DSA in most cosmic shocks. For example, even 
in the absence of any primordial seed field, magnetic fields are likely to be spontaneously 
generated at "external cosmic shocks" by some combination of the Biermann battery (Kul- 
srud et al. 1997; Ryu, Kang & Biermann 1998) and the Weibel instability (Schlickeiser 
& Shukla 2003). Whether such seed fields are initially strong enough to reduce diffusion 
times to cosmologically interesting values or not, quasi-linear plasma theory and simulations 
show that the streaming motion of suprathermal particles can induce resonant Alfven wave 
generation and strong field amplification upstream of collisionless shocks (e.g., Bell 1978; 
Quest 1988; Lucek and Bell 2000; Bell and Lucek 2001). Bell and Lucek (2001) have argued 
that magnetic fields can ultimately be amplified to produce magnetic pressures that are a 
substantial fraction of the ram pressure through the shock, PqU 2 s . That limit would corre- 
spond to field strengths ~ 4 x 10~ 7 /i(l + 5) l l 2 u s ^(\ + z) 3//2 Gauss, where h is the Hubble 
parameter in units of 100 km/sec/Mpc, 5 = Sp/p is the local overdensity ratio, z is the 
redshift of the epoch, and u s ^ is the shock speed in units of 10 3 km s -1 . A baryon density 
fraction, Q b = 0.04, is assumed. Even magnetic fields only one percent of this value would 
be sufficient to accelerate CRs to relativistic energies, and, hence establish nonlinear CR 
feedback on timescales of thousands of years in the shocks of interest. 

Once shocks develop nonlinear properties in response to CR feedback, there are several 
important characteristics that distinguish them from more familiar gasdynamic shocks: 1) 
CR shocks continue to evolve over relatively long times and broaden as they do so. While 
full thermalization takes place instantaneously at a simple, discontinuous jump in an "ideal" 
gasdynamic shock, CR acceleration and the corresponding modifications to the underlying 
flow depend on suprathermal particles passing back and forth diffusively across the shock 
structure. These processes develop, therefore, on the diffusion time scale, td, and diffusion 
length scale, Xd = tdU s , as alluded to earlier. Generally, these scales are expected to be 
increasing functions of particle momentum, so CR acceleration and shock evolution rates 
slow over time. For very high energy particles DSA time and length scales often approach 
the age and the radius of curvature of astrophysical shocks. 2) CR diffusion upstream of the 
shock discontinuity leads to strong pressure gradients in a shock precursor, enhancing the 
total compression through the shock transition over that in the "viscous" subshock. The 
total compression through a high-Mach-number CR shock can greatly exceed the canonical 
value r = (7^ + 1)/(7 S — 1) for strong gasdynamical shocks (e.g., Berezhko and Ellison 1999; 
Malkov, Diamond & Volk 2000) (where 7 9 is the gas adiabatic index). 3) Both the effective 
compressibility of the combined thermal-CR plasmas and the mean CR diffusion coefficient 
increase over time as relativistic CRs absorb more energy and as escaping CRs remove 
energy from the shock structure. Details of these properties depend on the CR momentum 
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distribution. Thus, downstream states of the CR modified shock cannot be found from 
simple "shock jump conditions", but rather have to be integrated time-dependently from 
given initial states or, for steady solutions, found in terms of some predefined limits to the 
CR spectrum (e.g., an upper momentum cutoff). 4) Energy transfer to the CRs rather than 
to the thermalized gas in the shock transition reduces the downstream thermal energy of a 
CR modified shock compared to a gasdynamic shock with the same shock speed and Mach 
number. The Mach number of the gas subshock is also reduced. Even in a shock with a very 
large total Mach number, the subshock Mach number and compression ratio are typically 
~ 3 (e.g., Malkov 1998; Berezhko and Ellison 1999; Kang, Jones, and Gieseler 2002). 

In order to study these processes closely we recently developed a numerical scheme that 
self-consistently incorporates a thermal leakage CR injection model based on the analytic, 
nonlinear calculations of Malkov (1998) (Gieseler et al. 2000) and implemented it into a 
combined Adaptive Mesh Refinement (AMR) gas dynamics and CR diffusion-convection 
code (Kang, Jones, and Gieseler 2002). As described in the following section, there is only 
one rather well constrained parameter needed to express this injection model; namely, the 
strength of nonlinear MHD waves immediately downstream of the plasma subshock with 
respect to the upstream longitudinal magnetic field strength. 

Previously, we applied our new code to a preliminary investigation of cosmic shocks 
emerging in large scale structure formation of the universe (Kang and Jones 2002; Kang 
2003). Utilizing the Bohm diffusion model in ID quasi-parallel, plane-parallel shocks we 
showed that the resulting CR modified shocks evolve roughly "self-similarly" if the shock 
structure is expressed in terms of diffusion time and length scales, and that the time asymp- 
totic states depend mainly on the shock Mach number, with only a weak dependence on 
the parameter needed in the injection model. We found for strong shocks that about 10~ 3 
of incoming thermal particles are injected into the CR population with our thermal leakage 
injection model over the long term, and up to ~ 1/2 of the shock kinetic energy flux mea- 
sured in the initial shock frame is transferred to CRs. In the present contribution we confirm 
and extend our previous studies by including simulations of shock models with a broader 
range of physical parameters suitable for application to cosmic structure formation shocks, 
and also by including the influence of a pre-existing, upstream CR population, since many 
of the shocks of interest will probably form in previously shocked plasma. Our numerical 
simulations here are limited to quasi-parallel shocks in which the mean magnetic field lines 
are parallel to the direction of shock propagation. In quasi-planar shocks this is not a strong 
limitation so long as the obliquity of the magnetic field is small enough to keep its intersec- 
tion subluminal and so long as the magnetic field is too weak to be dynamically important 
(e.g., Jokipii 1982; Drury 1983). In addition, since streaming motions inside cosmic sheets 
and filaments appear to stretch the field lines (e.g., Ryu, Kang & Biermann 1998), many 
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structure formation "internal" shocks are likely to be quasi-parallel. The principle concern 
may be the influence of obliquity on the injection of low energy CRs at the shocks, as we 
discuss below. 

In the following section we briefly outline our numerical methods in the CR/AMR 
hydrodynamics code. The simulation results are presented and discussed in §3, followed by 
a summary in §4. 



2. Numerical Calculations 

2.1. The CRASH code 

Our conservative, Eulerian CR/AMR hydrodynamics code, CRASH (Cosmic- Ray Amr 
SHock), solves the gasdynamic equations with CR pressure terms added for one dimensional 
plane-parallel geometry, along with the diffusion-convection equation for the CR momentum 
distribution function, g(p) = f(p)p 4 . The full CR shock transition includes a very wide 
range of length scales associated with the particle diffusion lengths that must be resolved 
to follow the shock evolution properly. To include CRs injected by thermal leakage, these 
scales must extend down close to the subshock thickness. In order to handle this dynamic 
range efficiently, an adaptive mesh refinement (AMR) technique is applied to regions around 
gasdynamic shocks, which are tracked as discontinuous jumps by a front-tracking method 
(Kang et al. 2001). An additional equation for the "Modified Entropy", S = Pg/p"' 9 ^ 1 , 
is solved to follow accurately the adiabatic changes outside of shocks, particularly in the 
precursor region of strong shocks (Kang, Jones, and Gieseler 2002). Readers are referred to 
these previous papers for further numerical details. 



2.2. The Bohm Diffusion Model 

DSA, along with the time and length scales for CR acceleration, and subsequent shock 
modification, depend on the spatial diffusion coefficient for the CRs. The diffusion coefficient 
can be expressed in terms of a mean scattering length, A, as n(x,p) = |Af, where v is the 
particle velocity. As noted earlier, for Alfven wave scattering one can write A ~ r g (B/5B) 2 
(e.g., Skilling 1975). The Bohm diffusion model, representing a saturated wave spectrum 
(SB ~ B) and the minimum diffusion coefficient, gives Kb = (l/3)r g v. If we assume "Bohm- 
like" diffusion with A = ( x r g with ( > 1, the diffusion coefficient can be expressed as 

«(M=«n(^) (p2 + 1)1/2 > (1) 
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where hereafter momentum, p, is expressed in units of mc and K n = ( x (l/3)mc 2 /(eB), 
represents the diffusion coefficient far upstream of the shock for CRs with p ^ 1.3. The 
assumed density dependence for k accounts for compression of the perpendicular component 
of the wave magnetic field and also inhibits the acoustic instability in the precursor of highly 
modified CR shocks (Drury and Falle 1986; Kang, Jones and Ryu 1992). Also, hereafter, we 
use the subscripts '0', '1', and '2' to denote conditions far upstream of the shock, immediately 
upstream of the gas subshock and immediately downstream of the subshock, respectively. 
Thus, p represents the far-upstream gas density. 

The so-called DSA diffusion time, td = k/u 2 s) defined in the introduction represents the 
time scale for hydrodynamical advection across the CR diffusion length, Xd, which, in turn 
represents the length on which CR diffusion upstream is balanced by advection towards the 
subshock. Conveniently, t d also provides a natural time unit to measure the acceleration 
time for individual CRs. The residence time for individual particles on each side of the 
shock is proportional to Xd/v, while the CR fractional momentum gain per shock crossing 
pair is proportional to Au/v. The resulting acceleration time scale for a particle to reach 
momentum p (e.g., Drury 1983) is defined as 

3 .Ki Ko, 8M 2 
Tacc(p) = - + - ~ TfT^tdip) ■ 2 

The approximate, Mach-number-based expression in equation (2) is based on the compression 
through a gasdynamic shock for 7 9 = 5/3 and sets u s = u±. Once again, td{p) = k{j>)/u 2 s . 
Thus, the CR acceleration time depends on both the speed and the Mach number (actually 
compression) of the shock in addition to the diffusion coefficient. In the limit of strong 
shocks (M s ^> 1), however, it becomes independent of the shock Mach number, and is 
about an order of magnitude greater than the nominal diffusion time, i.e., r acc pa 8t d . We 
note, however, that in CR modified shocks the compression factor across the total shock 
structure can be greater than that of strong gasdynamic shocks, and that higher energy 
CRs with greater diffusion length see a greater compression ratio. Consequently, the mean 
acceleration timescale at a fixed CR energy is somewhat smaller than that given in equation 
(2). 



2.3. The Thermal Leakage Injection Model 

In the "thermal leakage" model for CR injection at shocks, most of the downstream 
thermal protons are locally confined by nonlinear MHD waves and only particles well into 
the tail of the postshock Maxwellian distribution can leak upstream across the subshock 
(e.g., Ellison and Eichler 1985; Malkov 1998; Malkov and Drury 2001). In particular, "leak- 
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ing" particles not only must have velocities large enough to swim against the downstream 
flow in order to return across the shock, they must avoid being scattered by the MHD waves 
that mediate the plasma subshock. The latter effect further strongly filters particles against 
leaking. Thus, the ratio of the breadth of the postshock thermal velocity distribution to 
the downstream flow velocity in the subshock rest-frame is central to the injection prob- 
lem. Malkov and Volk (1998) and (Malkov 1998) developed a nonlinear analytic model to 
incorporate these features, and we have adapted that model to our simulations. In order to 
model this injection process numerically Gieseler et al. (2000) constructed a "transparency 
function", T esc (e,v) that expresses the probability of supra-thermal particles at a given ve- 
locity, v, leaking upstream through the postshock MHD waves. One free parameter, defined 
by Malkov and Volk (1998), controls this function; namely, e = Bq/B±, which is the ratio 
of the amplitude of the postshock MHD wave turbulence to the general magnetic field 
aligned with the shock normal, B . As noted by Malkov and Volk (1998) this parameter is 
rather well constrained, since 0.3 < e < 0.4 is indicated for strong parallel shocks. However, 
such large values of e lead to very efficient initial injection and most of the shock energy is 
quickly transferred to the CR component for strong shocks of M s > 30 (Kang and Jones 
2002), causing a numerical problem at the very early stage of simulations. Consequently, 
we reduced the shock transparency by setting e = 0.2 in this study. Kang (2003) showed 
that time asymptotic states of the CR shocks depend only weakly on the values of e in 
this range. On the other hand, smaller e values may also be more in tune with behaviors 
in moderately oblique shocks, since finite magnetic field obliquity also reduces shock trans- 
parency and the rate of thermal injection (e.g., Ellison, Baring & Jones 1995; Volk, Berezhko 
& Ksenofontov 2003). We note that Volk, Berezhko & Ksenofontov (2003) estimate that 
the injection rate for angles > 35° falls below the critical injection rate for efficient shock 
modification, £ ~ 6 x 10~ 5 (as defined in equation 5), and CR acceleration becomes very 
inefficient. 



2.4. Normalization of Physical Variables 

The ideal gasdynamic equations in ID planar geometry do not contain any intrinsic time 
and length scales, but in CR modified shocks the CR acceleration and the precursor growth 
can be characterized by the diffusion scales, td{p) and Xd{p). K n = k{j> 1.3) provides a 
useful canonical value, since nonlinear feedback from CRs to the underlying flow becomes 
significant typically by the time transrelativistic CRs are accelerated. Consequently, we 
normalize our shock evolution times and structure scales by t n = n n /u\, and x n = K n /u n , 
respectively, where u n is a characteristic flow speed. One expects intrinsic similarities in the 
dynamic evolution and structure of two CR shock models with the same Mach number, but 
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with different shock speeds, so long as the results are expressed in terms oft n and x n . This 
approach removes the need for any particular choice in n n , but we note for convenience that 
10 23 C/B _ 7 cm s 2 , where B = B_ 7 x 10 7 Gauss. Then for example, a characteristic 
diffusion time, t n = n n ju\ = 3.1 x 10 7 (/(B^ 7 u^ 3 ) sec, and x n = 3.1 x 10 15 (/(B_ 7 u nj3 ) cm, 
where u n = u Qt 3 x 10 3 km/sec. 

In the simulations described below we express pressures and energy densities in units 
P-n — Pn^ 2 , where p n will be an appropriate fiducial gas density. Again, for reference, we can 
express these quantities in practical units in terms of the preshock cosmic overdensity ratio, S, 
as p n w 7.5xl0~ 31 (l + 5)/i 2 (l + ,2) 3 gcm- 3 and P n « 7.5 x 10~ 15 /i 2 (l + 5)(l + ^) 3 M 2 i3 dyne cm" 2 , 
with the same cosmological choices mentioned in the introduction. 

When the shock speeds are less than a percent or so of the speed of light, freshly injected 
CRs are essentially nonrelativistic (p in j oc \f T P 2 oc (u s /c)), so both the diffusion coefficient and 
CR partial pressure take power law forms; thus, allowing self-similar behaviors in modified 
shock evolution. On the other hand, once the shock speed exceeds a percent or so of the 
speed of light, the full expression for the particle speed, v = p/ y 'p 2 + 1, must be applied, 
so neither the diffusion coefficient nor the CR partial pressure permit self-similar scalings. 
Consequently, the physical shock speed does influence the evolution of the fastest shocks we 
simulated in this study. 



2.5. Simulation Set Up and Model Parameters 

We model shocks that form in an accretion flow entering through the right boundary 
in a ID simulation box,[0, x max ]. The inflow has constant density, speed, and pressure. This 
flow reflects off the left boundary, or piston, to initiate the shock, which evolves in response 
to CRs accelerated via DSA as it propagates to the right. With the normalization constants 
described above, po = 1, Uq — — 1, and P gi o — (V7s)^o~ 2 m c °de units, where M is the 
accretion flow Mach number. The gas adiabatic index, 7 g = 5/3. The shock normalization 
constants are set by two parameters, M and c S;0 , through 

u Q = \u \ = c S)0 M , (3) 

where c Sj o — 15 km s^ 1 (To/10 4 ) 1//2 is the upstream sound speed of the accretion flow. We 
assume the preshock gas is fully ionized with a mean molecular weight, 0.61. We do not 
include complexities due to interactions between neutral and ionized gas in our numerical 
treatment. We also ignore added energy transport processes such as radiative cooling, the 
diffuse radiation emitted by the shock-heated gas and thermal conduction. However, they 
should be included in future studies for more realistic model calculations. 
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For a gasdynamic shock with compression ratio, r, the speed of the shock reflected from 
the piston is u s = u n /(r — l) in the piston (i.e., simulation) frame and u' s = u n r/(r — 1) in the 
upstream flow frame. From this relation, the initial shock speed with respect to the upstream 
flow spans the range u' s0 = (4/3 — 3/2)w n for accretion Mach numbers, 2 < M < 100. 
However, as CR pressure builds in the shock precursor, the total compression through the 
shock increases, and the shock slows down. That is, the instantaneous shock speed, u' s (t), 
relative to the accretion flow decreases over time. In addition, the gas subshock weakens as 
the flow is decelerated and heated in the precursor by the CR pressure gradient. 

We ran most simulations until time tf/t n = 20 — 40, which is generally sufficient to 
reach an approximate dynamical steady state, as discussed in §3. Several much longer 
simulations based on a new "coarse grained momentum finite volume" code in development 
have also been carried out to confirm the time asymptotic behaviors, as discussed in §3.5. 
The simulated space is x/x n = [0,20] for tf/t a = 20 and x/x n = [0,40] for tf/t a = 40 with 
N = (1 — 2) x 10 4 spatial zones on the base grid. The simulations were carried out on a 
base grid with Axq = (2 — 4) x 10~ 3 using Z max = 7 additional grid levels with a refinement 
of two, so Ax j = Axq/2 7 at the finest grid level. The number of refined zones around the 
shock is N r f = 50 on the base grid and so there are 2N r f = 100 zones on each refined 
level. To avoid technical difficulties the AMR technique is turned on and the CR injection 
and acceleration are activated only after the shock propagates to a half length of the pre- 
defined AMR refinement region in the base grid. This initial delay of the CR injection and 
acceleration should not affect the final outcomes. For all models we use 230 uniformly spaced 
logarithmic momentum zones in the interval \og(p/m p c) = [logp > l°gPi] — [—3.0, +3.0] 

Shocks formed by steady accretion flow onto cosmic structure pancakes provide a con- 
venient prototype of our modeled shocks. We, therefore, refer to the upstream flows as "ac- 
cretion flows" , and specify the flow speed and Mach number with respect to the reflecting 
symmetry plane. The models are also applicable to other cosmic structure shocks associated 
with filaments, knots and individual clusters. In Kang and Jones (2002) we considered a 
set of models with fixed infall speed, uq = —1500 km s~ x ("constant Mo" models), while the 
preshock temperature was varied according to T = 10 4 K(100/M ) 2 with 2 < M < 100. 
In Kang (2003), on the other hand, we set the preshock temperature at T = 10 4 K ("con- 
stant T " models), and varied the shock speeds according to u = (—15 km s _1 )M , with 
5 < M < 50. We now extend the study of Kang (2003) to include higher preshock tempera- 
tures, and also to include the effects of pre-existing CRs in the shock evolution. In particular 
we include: 1) T = 10 6 K shocks without pre-existing CRs, 2) T = 10 4 K shocks with pre- 
existing CRs, and 3) T = 10 6 K shocks with pre-existing CRs. For T = 10 4 K upstream 
conditions models with 5 < Mq < 5 ( 75 km s _1 < u n < 750 km s _1 ) are considered, while for 
T = 10 6 K conditions we include cases with 2 < M < 50 (300 km s _1 < u n < 7500 km s _1 ). 
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We note, for T = 10 6 K conditions, models with M < 20 (u n < 3000 km s^ 1 ) would be 
relevant for cosmic structure shocks, but models with M = 30 and 50 are included to explore 
the similarity issue discussed above in §2.4. In determining these sets of model parameters, 
we considered the following factors. In astrophysical environments photoionized gas of 10 4 
K is quite common. For example, Ly a clouds and the intergalactic medium in cosmological 
sheets and voids are believed to be photoionized during the reionization epoch (Loeb and 
Barkana 2002). Also the preshock gas can be photoionized and heated to ~ 10 4 K by the 
diffuse radiation emitted by hot postshock gas when u s > 110 km s -1 (Shull and McKee 
1979; Kang and Shapiro 1992; Draine and McKee 1993). On the other hand, hot and ionized 
gas of > 10 6 K is also found in the hot phase of the ISM (McKee and Ostriker 1977; Spitzer 
1990) and in the intracluster medium of X-ray clusters (Fabian 1994). 

Since cosmological simulations indicate that most diffuse plasma in large scale cosmic 
structures is likely to be shocked more than once during structure formation, it is likely 
that many "internal shocks" process plasma that carries previously accelerated CRs. We, 
therefore, include shock simulations including these populations. For pre-existing, upstream 
CRs we assumed a simple power-law spectrum smoothly connected to the upstream thermal 
Maxwellian distribution (/m) at a momentum p M as 

fu P (p) = !m (pm ) (p/pm ) ~ qin (4) 

for p M < p < Pi, where pu = 5(m p /csTo) 1//2 , p\ = 10 3 . For T = 10 4 K shocks we considered 
in the pre-existing CR cases an upstream CR pressure condition, P Cj0 = 0.25P gi o, using 
q in = 4.7. For T = 10 6 K shocks we included three different upstream conditions; namely, 
P c , = 0.25P gi0 (<?;„ = 4.5), P Ci0 = 0.5P g , (<? in = 4.4) and P c , = P g , (q in = 4.3). The CR 
spectral forms represent test-particle populations for previous shocks with moderate Mach 
numbers in the range, 2.5 < M s < 4. Table 1 provides a summary of the new models 
computed for the present study. 



2.6. Injection and Acceleration Efficiencies 

It is useful to characterize the efficiency with which low energy CRs are injected at the 
shocks. We define the injection efficiency as the fraction of particles that have entered the 
shock from far upstream and then injected into the CR distribution: 

fX 



Jln u' s (t')dt' 



m = - — T :, : (5) 



where /cr is the CR distribution function, n Q is the particle number density far upstream, 
and U is the time when the CR injection/acceleration is turned on (t, ss 3t n for M > 5 and 
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U « 2t n for M = 2). 

If the subshock speed becomes steady, then n u' s is the same as the particle flux swept 
by the subshock, niW su b, where n\ is the particle number density immediately upstream of 
the subshock. However, in our simulations these two flux measures can differ by up to 10% 
for strong shock models, because the shock speed is changing slowly. 

As a measure of efficiency of CR energy extraction at shocks, we define the "CR energy 
ratio", $; namely, the ratio of the total CR energy within the simulation box to the kinetic 
energy in the initial shock frame that has entered the simulation box from far upstream, 

= r"d* gcR( M) 

1 ' 0.5p„« )3 t 1 1 

All shock models have the same normalized upstream mass density and velocity, po and 
u , but different upstream gas thermal energy density, -P g ,o/(7 g — 1), depending on M . 
Therefore, we use the kinetic energy flux rather than the total energy flux to normalize 
the "CR energy ratio". Alternatively, the ratio defined with the total energy flux, F tot = 
w' s [0.5po(X,o) 2 + IgPgfl/ilg - 1) + 7c,o-P c ,o/(7c,o - 1)], can be calculated from the following 
relation, 

^tot{t) « — 3 (7) 

where we assume 7 Cj0 ~ 5/3 for the preshock CR population. This shows that $ to< ks $ for 
large values of M , while § tot « (3/5)$ for M = 2 and P Cj0 = P gfi . 



3. Simulation Results 



We can define three key stages of shock evolution in response to CR feedback: 1) De- 
velopment of the shock precursor that slows and heats the flow entering the gas subshock, 
reducing the Mach number of the latter; 2) achievement of an approximately time-asymptotic 
dynamical shock transition, including nearly steady postshock CR and gas pressures; 3) con- 
tinued, approximately "self-similar" evolution of the shock structure for Bohm-type diffusion, 
as CRs are accelerated to ever higher momenta. Once stage (3) is reached, there is little 
change in the dynamical properties of the shocks, and, in particular, little change in the 
total efficiency with which kinetic energy is transferred to CRs. Thus, given the rapidly 
increasing computational cost to simulate acceleration of CRs to still higher energies with 
Bohm diffusion (which more than quadruples for each factor of two increase in maximum 
CR momentum), we deemed it sufficient to the questions being addressed to terminate the 
simulations at the times indicated. 
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As Kang (2003) showed, the evolution of CR modified shocks with Bohm-type diffusion 
are mainly determined by the shock Mach number for a given injection parameter. This 
includes the efficiency of energy transfer to CRs, $. That is, two shocks with the same 
Mach number, but with different shock speeds evolve qualitatively similarly when expressed 
in terms of characteristic diffusion scales, t n and x n . During the early evolutionary stage, 
however, the similarity between the two shocks may break down as pointed out in §2.4. 
As described there, this happens if the shock speed is larger than a percent or so of the 
speed of light, on account of the different momentum- velocity scalings for nonrelativistic 
and relativistic particles at injection. As a result, the ratio of P c /(p n u^) increases faster and 
nonlinear modification sets in earlier in terms of t/t n in the shock with lower speed. Thus, 
depending on when the transition from nonrelativistic to relativistic population occurs, early 
evolution of the two shocks can be different (e.g., see M = 30 models in figure 2), although 
they still eventually approach similar time-asymptotic states. 

The other factor that can affect the similarity is the presence of pre-existing CRs. For 
strong shocks of M s > 10, the preshock CR pressure is only a few percent of the total energy 
flux even when P c ~ P g far upstream. In those cases the presence of pre-existing CRs is 
effectively equivalent to having a slightly higher injection rate (i.e., larger e), which speeds 
up the initial shock evolution. For strong shocks, our previous studies have shown that the 
time asymptotic CR acceleration efficiency depends only weakly on the injection rate (Kang, 
Jones, and Gieseler 2002; Kang and Jones 2002), provided the injection remains above a 
minimum threshold for "efficient" DSA (see also, e.g., Berezhko et al. 1996). On the other 
hand, for weak shocks, a pre-existing CR pressure comparable to the upstream gas pressure 
represents a significant fraction of the total energy entering the shock, so pre-existing CRs 
have far more impact. Also the time asymptotic CR acceleration efficiency in weak shocks 
depends sensitively on the injection rate and increases with e. Hence we should expect that 
relatively weak CR shocks (M s < 5) will be substantially altered by the presence of a finite 
upstream P c . 



3.1. Evolution of Shock Structures 

Figure 1 compares the time evolution of two CR modified shock models with accretion 
Mach number M = 5, but with different physical flow speeds and upstream CR conditions. 
The light lines show the flow model with T = 10 4 K and the pre-existing CR population, 
f up oc p -4 ' 5 , while the heavy lines show the model with To = 10 6 K and f up oc p~ 4 ' 7 . The 
pressure due to the pre-existing CR particles is P Ci0 = 0.25P 9j0 in both models. Note that the 
physical time and length scales are different, since t n oc u~ 2 and x Q oc u' 1 (u n = 75 km s" 1 
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for T = 10 4 K models and u n = 750 km s^ 1 for T = 10 6 K models). 

In the model with u D = 75 km s _1 CR particles are injected at much lower injection 
momenta. Since for a fixed Mach number their initial acceleration times scale as r acc /t n oc 
p? n j, the acceleration and shock modification takes place more rapidly in these time units 
compared to the model with u n = 750 km s _1 . The value of P c /(p I1 u 2 l ) approaches 0.8 for 
the model with u n = 75 km s -1 , while it approaches 0.65 for the model with u n = 750 km s _1 . 
Otherwise, the overall evolution is roughly similar. 

Figure 2 shows that even for much higher Mach numbers, M = 30, two shocks with u n = 
450 km s _1 and u n = 4500 km s -1 approach similar time-asymptotic states, although their 
early evolutions are very different. However, as the shock speed increases further, substantial 
differences develop in the early evolution. For M = 50 models, for example, unlike the 
shocks shown in the figures, the overall shock structures then evolve much more slowly in 
our normalized time units, since the injection momentum is no longer nonrelativistic. Then 
the Bohm-like diffusion coefficient and the CR partial pressure near this momentum do not 
scale with p 2 as they do for nonrelativistic momenta. The approximate self-similarity of 
the shock evolution breaks down under those conditions. Even those shocks, however, do 
approach similar time asymptotic states for a given Mach number. 

The flow entering the piston is decelerated first by the precursor and then by the sub- 
shock in our simulation. For a pure gas fluid with a constant adiabatic index, the postshock 
flow speed relative to the piston should be uniformly zero. This should be approximately 
true for the models with inefficient CR acceleration in our simulations (see, for example, 
Fig. 4 below). In CR modified shocks, however, the "effective" adiabatic index decreases due 
to increasing CR energy, so the total compression through the shock increases and the shock 
slows down over time. This transition in the compressibility of the flow leads to a net forward 
flow velocity of 0.04 < u/u n < 0.07 in the piston rest frame for the two models shown in 
Figs. 1 and 2. However, the net flow speed is decelerated by the positive gradient of the 
combined postshock pressure (P c + P g ) and approaches a smaller value (0.02 < u/u n < 0.03) 
at the terminal simulation time. 

As noted in previous time- dependent studies (see the references in §1), the compression 
factor for these CR dominated shocks relaxes to values still higher than that for a gasdynamic 
shock, even a for a fully relativistic gas, even though energy is not leaving the system in our 
simulations. The reason has to do with the evolving distribution of energy in the system. In 
particular, for a diffusion coefficient increasing with p, especially for high Mach number, CR 
dominated shocks, the compression factor can remain much higher than what is expected in 
a steady-state, energy-conserving shock jump for relativistic gas (i.e., p2/rhoo = 7). This 
is because energy is being continuously transported out of the region of the subshock by 
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diffusion of the highest energy CRs with an ever-increasing diffusion scale Xd(p max ). More 
extensive discussion can be found in Kang, Jones, and Gieseler (2002). 

3.2. Effect of Pre-existing CRs for Weak Shocks 

Figure 3 compares the evolution of shocks with and without pre-existing, upstream CRs. 
Four sets of shocks are illustrated: 1) M = 5, u n — 75 km s -1 , and T = 10 4 K (upper 
left panel), 2) M = 5, u n = 750 km s _1 , and T = 10 6 K (upper right panel), 3) M = 30, 
u n = 450 km s" 1 , and T = 10 4 K (lower left panel), and 4) M = 30, u n = 4500 km s _1 , 
and T = 10 6 K (lower right panel). The approximate self-similarity for the evolution of 
these shocks, when described in normalized units is apparent. In addition one can see that 
the presence of pre-existing CRs is clearly more important for low Mach number shocks. 

It is especially important to consider the influence of pre-existing CRs in low Mach 
number shocks, since internal shocks with M s ~ 3 may be the most important for dissipation 
in shock-heated gas inside nonlinear structures (Ryu et al. 2003). We already know that CR 
acceleration efficiency in low to moderate strength shocks depends sensitively on the injection 
rate and increases with e (Kang and Jones 2002). Our new simulations show that the CR 
acceleration efficiency in these shocks depends similarly on the pre-existing (upstream) CR 
pressure. That trend was apparent in figure 1. Figure 4 emphasizes the point impressively 
for even lower Mach number shocks. It shows the time evolution of the models with accretion 
flow Mach number Mq = 2 and uq = —300 km s _1 , which drive shocks with initial Mach 
number of M s = 3 into the preshock gas with To = 10 6 K. Three models are compared in 
the figure: 1) a model with no pre-existing CRs (solid lines), 2) a model with P c0 = 0.25P 9i o 
(dotted lines), and 3) a model with P c0 = 0.5P Si o (dashed line). For the model without the 
pre-existing CRs and with injection parameter, e = 0.2, the CR injection is small and the 
resulting CR acceleration is inefficient. The solution is essentially a test-particle solution, 
with minimal CR modifications to the flow structure. There are a couple of related reasons 
for the impact of the upstream CRs in the two other models. First, the pressure of the 
pre-existing CRs cannot be neglected in the shock transition compared to the total shock 
energy flux. In addition, the inflowing CRs act like a large source of injected CRs, effectively 
enhancing the injection rate (see also Blasi 2003). Kang and Jones (2002) also considered 
similar M = 2 models without the pre-existing CRs (u = —1500 km s _1 , and T = 2.5 x 10 7 
K), but with a range of injection parameters, e = 0.2 — 0.4. We found that the postshock 
CR pressure increased over this range of e from P c = 0.04 to P c = 0.24, while the associated 
injection rate rose from £ = 10 -4 to £ = 10~ 2 (see figure 8 in Kang and Jones 2002). These 
results imply that even in low Mach number shocks the CR acceleration can have significant 
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dynamical influences, when the injection rate is as high as £ ~ 10~ 2 or when the upstream 
CR pressure is greater than a few times 10 % of the gas pressure. This could be important in 
the context of cosmic structure shocks, especially for "internal shocks" in which the upstream 
plasma already carried a significant population of CRs produced in an earlier era. 

3.3. The Particle Distribution 

Time evolution of the distribution function g{p) = f{p)p A at the gas subshock is shown 
in figure 5 for the models with To = 10 6 K and with f up oc p~ 4 ' 5 . The dotted line shows 
the initial Maxwellian distribution with T and the specified power-law distribution for the 
pre-existing CR population. The postshock thermal distribution during the early evolution 
(before shock modification) should have Ti oc M5 2 , so the injection momentum, p-^ oc Mq 
for t/t-a < 1. However, as energy transfer to CRs increases over time and the gas subshock 
weakens, the peak in the Maxwellian distribution shifts to lower momenta. The model with 
M = 30 shows this effect clearly. This nonlinear feedback is greater for higher Mach number, 
so that the postshock temperature of evolved CR modified shocks depends on M much more 
weakly than the standard M5 2 . 

It is useful to identify p max as the momentum above which g{p) drops sharply, character- 
izing the effective upper cutoff in the CR distribution, at a given time since shock formation. 
This momentum is generally easily identified. We find that it is approximately related to 
the age of the shock as p max ~ 0.4(t/t n ) for p max > 1. This explains why values of j9 max are 
similar at a given value of t/t Q for all models considered here, independent of u n , injection 
momentum, or shock Mach number. The previously described back reaction of the CRs on 
the gas subshock acts to move the shock towards a limiting form and postshock state. On 
the one hand, a fraction of CR particles injected at early times continues being accelerated 
to higher momenta, so that p max (t) always increases until they are subject to escape due 
to geometry or insufficient scattering. On the other hand, particle injection becomes less 
efficient as the subshock feature weakens. Because of this reduced source spatial diffusion 
causes the number density of CRs near the subshock at a given momentum to decrease with 
time for any p < p max (t). Combined, these effects cause the postshock CR pressure, P c , 
to approach a time-asymptotic value, and, for Bohm-like diffusion, the shock structure to 
become approximately self-similar. 
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3.4. Cosmic Ray Injection and Acceleration Efficiencies 

Figure 6 illustrates the time averaged injection efficiency, £(£), defined by equation (9) for 
models without pre-existing CRs. The calculated injection efficiency ranges from ~ a few x 

to ~ 10 -3 , depending on the subshock Mach number, the shock speed, and the injection 
parameter, e. As outlined in §2.4 and discussed fully in Kang, Jones, and Gieseler (2002), 
the modeled thermal-leakage injection process is less efficient for weaker subshocks and for 
smaller e (stronger wave trapping of suprathermal particles in the postshock environment). 
Injection is also less efficient for lower shock speeds, because the diffusive flux of injected 
particles crossing the subshock is smaller for smaller injection momentum (p in j oc u s ). So, 
for a given Mach number the models with higher shock speed (i.e., higher T in figure 6) 
have larger values of £(£). For strong shocks the subshock weakens significantly due to the 
CR nonlinear feedback, so the injection efficiency slowly decreases over time. 

Figure 7 shows for the models with T = 10 4 K the CR energy ratio, $, the CR pres- 
sure and the gas pressure at the subshock normalized by the far upstream ram pressure 
in the initial shock frame, i.e., P c ^/ (pqu 1 ^, and P g ,2/(Po M s o)- I n addition the bottom 
panels show the CR pressure normalized by the far upstream ram pressure in the instanta- 
neous shock frame, P c ^/ (pou' s 2 ), The left panels illustrate the models without pre-existing 
CRs, while right panels show the models with P Cj0 = 0.25P 9i0 - As described in the previ- 
ous subsection, the postshock P c for all Mach numbers increases until a balance between 
injection/acceleration and advection/diffusion of CRs is achieved. After that time, which 
corresponds in these simulations to the generation of moderately relativistic CRs in the 
shock, the postshock CR pressure remains almost steady. Simultaneously, the CR energy 
ratio, $, also initially increases with time, but then asymptotes to a constant value, once 
P Cj2 has reached a quasi-steady value. The asymptotic values result from the "self-similar" 
behavior of the P c distribution. Since both the numerator and denominator in equation (6) 
increase approximately linearly with time in the case of the self-similar evolution of E c , $ 
becomes steady when P c2 becomes constant. As illustrated in Figure 7, time- asymptotic 
values of $ increase with the accretion Mach number M , approaching $ xs 0.5 for M > 20 
by the simulation termination times. As expected from the discussion in §3.4, the presence of 
pre-existing CRs increases the asymptotic postshock P c for the weaker shock, Mo = 5 model, 
while it does not affect the strong shock models with Mq > 10. Figure 7 also illustrates the 
well-known result that strong shocks can be mediated mostly by CRs while the postshock 
gas thermal energy is reduced greatly compared to that in conventional gasdynamic shocks 
of the same Mach number. As a result of nonlinear feedback the shock slows down and so 
the ratio of P c ,2/(pou' s 2 ) can reach ~ 0.9, while the ratio of P Cj2 / '(pou' s ^) reaches ~ 0.6 in 
these models. 
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In figure 8 we present the same quantities for the models with T = 10 6 K. As explained 
in the previous section, due to the higher injection momenta, the CR acceleration is slower 
to develop for the models with M > 20, compared to the same Mach number model with 
lower shock speed. Especially for the M = 50 model, the postshock P c continues to increase 
at the termination time (t/t Q = 40) of the simulation. Again, the presence of pre-existing 
CRs significantly increases the postshock P c at low Mach number models (M < 5). We 
show the same quantities for low Mach number simulations with even higher pre-existing P c 
in figure 9. We see that the postshock P c increases with the pre-existing P c for these low 
Mach number models. 



3.5. Time Asymptotic Behaviors 

As discussed above, some characteristics of the CR modified shocks seem to have reached 
quasi-steady states well before the termination time for most of the models considered in 
this work. We have further explored these time asymptotic behaviors by performing much 
longer simulations with a new "Coarse-Grained Momentum Finite Volume" (CGM) numer- 
ical scheme under development. CGM utilizes the fact that the distribution function can 
be approximated as a piecewise power-law within a given broad momentum bin and so only 
roughly one momentum bin per octave is necessary to follow the CR population numerically. 
Details and tests of the CGM numerical method will be presented in a forthcoming paper 
(Jones and Kang 2004). In figures 10 and 11, we show the simulation results of Mo = 5 and 
M = 30 models, respectively, which were calculated by this CGM method up to t/t n = 1000. 
No pre-existing CRs were assumed and 14 momentum bins were used for these calculations. 
The heavy solid lines are for the simulation results at t/t n = 20 for the M = 5 model and 
t/t n = 40 for the M = 30 model which were calculated by the "fine-grained momentum" 
method described in §2 and used for all the other results presented in this paper. Those 
simulations were carried out with 240 momentum bins. These figures confirm very firmly the 
aforementioned "self-similarity" of the shock evolution and the time-asymptotic behaviors. 

Thus we show in figure 12 values of the CR energy ratio $ at the termination time as 
a function of the initial unmodified shock Mach number, M s , as a simple measure of the 
CR acceleration efficiency For convenience in comparing with other plots Values of M s are 
listed against M in Table 1. Flows with T = 10 4 are shown in the upper panel, while 
shocks developed in T = 10 6 flows are represented in the lower panel. For comparison, we 
also plotted in the upper panel results from "constant uq" models with uq = —1500 km s _1 
and T = 10 4 (100/M ) 2 K, originally presented in figure 7 of Kang (2003). Those models, 
which did not include a pre-existing CR population, behave similarly to the new models 
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considered in the present work. In each series of simulations the efficiency asymptotes to 
$ 0.5 for large M s . This asymptotic efficiency is independent of P Cj o/-P 9 ,o; i-e., the level of 
pre-existing CRs. It also reflects the fact that the postshock CR pressures in the high Mach 
number shocks evolve in all cases roughly to 60% of the inflowing ram pressure measured in 
the initial shock rest frame, poit' s q. In weak to moderate strength shocks, on the other hand, 
the CR efficiencies shown in figure 12, or the postshock CR pressures illustrated in figures 7 
and 8, depend strongly on P c ,o/P g ,o- A rough dividing line for these two behavior domains 
is M s ~ 20. For the M s = 3 model without the pre-existing CRs and e = 0.2, the injection 
rate ^«2x 10~ 5 and <3> approaches to only 0.04. Since 0.3 < e < 0.4 is indicated for strong 
shocks (Malkov and Volk 1998) and weaker wave fields are expected in lower Mach shocks, 
leading to larger values of e, we show three addition models for e = 0.25, 0.3, and 0.4 with 
M s = 3, To = 2.5 x 10 7 K and no pre-existing CRs, originally presented in figure 8 of Kang 
and Jones (2002), in the upper panel of figure 12. In Ryu et al. (2003) the CR efficiency 
model, $(M S ), for e = 0.3 was adopted in the calculation of the CR energy generated by 
cosmic stricture shocks. We note that the functional minimum in $(M S ) for the models with 
pre-existing CRs is an artifact of the definition of $, which is normalized by the inflowing 
kinetic energy, not the total energy. The thermal energy flux is a significant fraction of the 
kinetic energy flux for low Mach number shocks. There is no minimum in the function $ tot , 
as defined in equation 7. For instance the value $(M S = 3) ^ 0.67 shown in figure 12, 
becomes & to t(M s = 3) ^ 0.4. 

As also pointed out in §3.2, the increases in CR energy ratio for weak shocks when they 
process plasma that contains pre-existing CRs, could be important to understanding the 
properties of cosmic structure formation shocks. In addition, most baryonic matter in clusters 
has been shocked more than once before the current epoch (Ryu et al. 2003). This means that 
many of the shocks primarily responsible for energy dissipation during structure formation 
may be influenced by CRs accelerated in previous shock events, since magnetic field strengths 
at levels discussed in the introduction should tie CRs up to moderately relativistic energies 
to the thermal plasma over cosmic time scales (e.g., Volk, Aharonian & Breitschwerdt 1996). 
In that case, those shocks may be substantially more efficient at CR acceleration than one 
might otherwise expect. 



4. Summary 

Energy dissipation at collisionless shocks involves complex physical processes and plays 
crucial roles in the evolution of shock-heated astrophysical plasmas. Here our focus is CR 
acceleration via the first-order Fermi process at cosmic shocks that form in and around 
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nonlinear structures during the formation of large scale structure in the universe. Using 
our cosmic-ray AMR Shock code (Kang, Jones, and Gieseler 2002), we have calculated 
the CR acceleration at ID quasi-parallel shocks driven by plane-parallel infall flows with 
| Mo I — 75 — 7500kms _1 and preshock temperatures of T = 10 4 — 10 7 6 K. Mach numbers of 
the resulting shocks range over 2 < M s < 133, which covers the characteristic shock Mach 
numbers found in cosmological hydrodynamic simulations of a ACDM universe (Ryu et al. 
2003). 

The main purpose of this study is to explore how the CR acceleration efficiency at cosmic 
structure shocks depends on the preshock temperature, shock speed and pre-existing CRs as 
well as shock Mach number. This, in turn, affects the thermal and dynamical history of the 
shocked-heated baryon gas, since the energy transfer to CRs reduces the thermal energy of 
the gas and increases the compressibility of the gas flow, much like radiative shocks. 

Unlike pure gasdynamic shocks, downstream states of the CR modified shocks cannot 
be found from "shock jump conditions", so they have to be integrated time-dependently 
from given initial states or found from nonlinear analytical methods. Time-dependent simu- 
lations are computationally expensive, because the CR diffusion and acceleration processes 
contain a wide range of length and time scales. Fortunately, there are approximate similar- 
ity properties that can be employed to study the time asymptotic behavior of evolved CR 
modified shocks when they are mediated by Bohm-like diffusion. Firstly, the CR pressure 
approaches a steady-state value in a time scale comparable to the acceleration time scales 
for mildly relativistic protons after which the evolution of CR modified shocks becomes ap- 
proximately "self-similar" . This feature enables us to predict time asymptotic values of the 
CR acceleration efficiency, although we followed the CR acceleration up to only moderately 
relativistic energies {i.e., p max /(mc) ~ 16). We have confirmed these time asymptotic behav- 
iors through much longer simulations with a new, more efficient, numerical scheme (Jones 
and Kang 2004), which were extended to achieve the cut-off momentum p max /( m c) ~ 400. 
Secondly, two models with the same Mach number, but with different accretion speeds show 
qualitatively similar evolution and structure when the dynamical evolution is presented in 
terms of characteristic diffusion scales, t n = n n /u\ and x Q = K n /u n . Such similarities are only 
approximate, however, because the partial CR pressure and the Bohm diffusion coefficient 
for transrelativistic CRs do depend on the momentum m p c. Since the effective injection 
momentum is Pinj/tripC oc (u s /c), the initial evolution depends on the shock speed as well as 
Mach number. 

Finally, we state the main conclusions of our study, which can be summarized as follows: 

1) Suprathermal particles can be injected efficiently into the CR population at quasi- 
parallel cosmic shocks via the thermal leakage process. For a given injection parameter 
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defined in the text, e, the fraction of injected CRs increases with the subshock Mach number, 
but approaches £ ~ 10~ 3 in the strong shock limit. 

2) For a given value of e, the acceleration efficiency increases with the shock Mach 
number, but approaches a similar value in the strong shock limit. Time asymptotic values of 
the ratio of CR energy to inflowing kinetic energy, $, converge to $ « 0.5 for M s > 30 and 
it is relatively independent of other upstream or injection parameters (see figure 12). Thus, 
strong cosmic shocks can be mediated mostly by CRs and the gas thermal energy can be up 
to ~ 10 times smaller than that expected for gasdynamic shocks. 

3) For weak shocks, on the other hand, the acceleration efficiency increases with the 
injection rate (or e) and the pre-exiting P c (see figure 12). For cosmologically common 
M s = 3 shocks, for example, $ ~ 0.04 for e = 0.2 and $ ~ 0.2 for e = 0.3 in the absence of 
pre-existing CRs. The presence of a pre-existing CR population acts effectively as a higher 
injection rate than the thermal leakage alone, leading to greatly enhanced CR acceleration 
efficiency in low Mach number shocks. We found even for weak shocks of M s = 3, that up 
to 40 % of the total energy flux through the shocks can be transferred to CRs, when the 
upstream CR pressure is comparable to the gas pressure in the preshock flow. 
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M n = 5, e = 0.2, t/t = 2, 5, 10, 20 
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Fig. 1. — Time evolution of the shocks driven by ID accretion flows with Mq = u n /c Sy o = 5 
is shown at normalized times, t/t n = 2, 5, 10, and 20. The accretion flow enters from 
the right boundary (x/x n = 20) and is reflected at the piston (x/x n = 0), leading to a 
shock with M s = 6.8 propagating to the right. The flow speed is shown in the piston rest 
frame, so the flow is almost at rest near the piston. The leftmost profile corresponds to 
the earliest time. Light lines represent a flow with u n = 75 km s _1 , T = 10 4 K, and a 
pre-existing CR population with f up oc (p/pm)~ 4 ' 7 ■ Heavy lines represent the model with 
u n = 750 km s -1 , T = 10 6 K, and f up oc (p/pm)~ 4 ' 5 - For both models the upstream 
CR pressure is P c ,o/P g ,o = 0.25. The normalization diffusion time scale, t n = k d /m^, and 
diffusion length, x n = K n /w n are defined by the accretion speed of each model. The inverse 
wave-amplitude parameter e = 0.2 is adopted for both models. 
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Fig. 2. — Same as Figure 1 except M = 30. Light lines are for the model with u n = 
450 km s _1 , T = 10 4 K, and a pre-existing CR population of f up oc (p/^m) -4 ' 7 , while 
heavy lines represent the model with u n = 4500 km s _1 , T = 10 6 K, and f up oc (p/pm)~ 4 ' 5 - 
For both models the upstream CR pressure is P c ,o/Pg,o = 0.25. Note for T = 10 6 K model 
additional curves (heavy solid lines) are shown at t/t n = 40. 
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Fig. 3. — Time evolution of the CR pressure for the following models: 1) Mq = 5 and 
T = 10 4 K (upper left panel), 2) M = 5 and T = 10 6 K (upper right panel), 3) M = 30 
and T = 10 4 K (lower left panel), and 4) M = 30 and T = 10 6 K (lower right panel). 
Heavy lines are for the models without pre-existing CRs, while light lines are for the models 
with P Cj0 = 0.25P 5i0 . 
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M n = 2, £ = 0.2, t/t = 10, 20, 30, 40 
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Fig. 4. — Same as Figure 1 except M = 2, T = 10 6 K and u n = 300 km s _1 . The solid 
lines are for the model with P Cj0 = 0.0, the dotted lines for the model with P c = 0.25P ffi o, 
and the dashed lines for the model with P c0 = 0.5P 9)0 . 



-28- 




Fig. 5. — Evolution of the CR distribution function at the shock, represented as g = p 4 f(p), 
is plotted for the models of M = 5, 10, 30, and 50 with T = 10 6 K and f up oc Qo/pm) -4 ' 5 - 
The CR spectrum of the preshock flow is represented by the dotted line. For all models 
shown here the upstream CR pressure is P c ,o/P g ,o = 0.25. 
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Fig. 6. — Time averaged injection efficiency, £(t), for models without pre-existing CRs. Left 
panel shows the models with M = 5 — 50 and T = 10 4 K, while right panel shows the 
models with M = 2 - 50 and T = 10 6 K. 
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Fig. 7. — The ratio of total CR energy in the simulation box to the kinetic energy in the initial 
shock rest frame that has entered the simulation box from upstream, the postshock 

CR pressure, P c ,i y and gas pressure, P 9t 2, in units of upstream ram pressure in the initial 
shock frame, pou' s \ and the postshock CR pressure in units of upstream ram pressure in 
the instantaneous shock frame, Pqu' s 2 Left panels show the models with M = 5 — 50 and 
T = 10 4 K but without pre-existing CRs. Right panels for the same models but with the 
upstream CR pressure of P Cj o = 0.25P 9> o and f up oc (p/pju) -4 ' 7 - 
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Fig. 8. — Same as figure 9 except that the models with T = 10 6 K are shown. Left panels for 
the models without pre-existing CRs, while right panels for the models with P Cj0 = 0.25P 9i o 
and f up oc (jo/Pm)~ 4,5 • 
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Fig. 9. — Same as figure 9 except that the low Mach models with To = 10 6 K are shown. 
Left panels show the models with pre-existing CRs of P Cj0 = 0.5P g and f up oc (p/pm)~ 4A , 
while right panels for the models with P c0 = P g0 and f up oc (jo/pm) -4 ' 3 - Solid line is for 
Mo = 2 model, dotted line for M = 3 model, and dashed line for M = 5 models. 



-33- 



M =5 piston shock, k=k b , t= 20, 200, 400, 600, 800, 1000 
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Fig. 10. — Time evolution of the Mo = 5 model without pre-existing CRs calculated with 
"Coarse-Grained Momentum Volume" method is shown up to t/t n =1000. The leftmost 
profile corresponds to the earliest time, t/t n = 20. For comparison, the results of "Fine- 
Grained Momentum Volume" simulation are shown at the same time by the heavy solid line. 




Fig. 11. — Same as Figure 10 except M = 30. The results from two simulations are 
compared at t/t n = 40. 
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Fig. 12. — The CR energy ratio, $, at the simulation termination time of our simulations 
as a function of the shock Mach number, M s . Upper panel: models with T = 10 4 K, 
u s = (15 km s -1 )Ms, e = 0.2, and no pre-existing CRs (solid line with open circles), models 
with T = 10 4 K, u s = (15 km s~ l )M s , e = 0.2, and P cfl = 0.25P 9 , (solid line with filled 
triangles). For comparison we also show "constant uq models with \uq\ = 1500 km s _1 , 
T = 10 4 K(100/M ), e = 0.2, and no pre-existing CRs in the upper panel (dashed line with 
filled triangles). Three filled circles are labeled with the value of e = 0.25, 0.3, 0.4 and for 
the models with M s = 3, To = 5 x 10 5 K and no pre-existing CRs. Lower panel: models with 
T = 10 6 K, u s = (150 km s _1 )M s , e = 0.2, and a various level of pre-existing CR pressure. 
See Table 1 for the relation between the accretion Mach number M and M s . 
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Table 1. Model Parameters for Flows Upstream of Shocks 



T 


Kl 


M a 




Pc,0/Pg,0 


Qin 


(K) 


( km s _1 ) 










10 4 


(15)M 


5, 10, 20, 30, 50 


6.8, 13.4, 26.7, 40., 66.7 


0.0 




10 4 


(15)M 


5, 10, 20, 30, 50 


6.8, 13.4, 26.7, 40., 66.7 


0.25 


4.7 


10 6 


(150) M 


2, 3, 5, 10, 20, 30, 50 


3., 4.2, 6.8, 13.4, 26.7, 40., 66.7 


0.0 




10 6 


(150)M 


2, 3, 5, 10, 20, 30, 50 


3., 4.2, 6.8, 13.4, 26.7, 40., 66.7 


0.25 


4.5 


10 6 


(150)M 


2, 3, 5, 10 


3., 4.2, 6.8, 13.4 


0.5 


4.4 


10 6 


(150)M 


2, 3, 5, 10 


3., 4.2, 6.8, 13.4 


1.0 


4.3 


10 S M 2 


1500 


1.5, 2, 3, 5, 10, 30, 100 


2.4, 3., 4.2, 6.8, 13.4, 40., 133. 


0.0 





a Mach number of piston with respect to inflow. 

b Mach number of initial gas shock with respect to infall. 

c Power-law index of upstream CR population, i.e., f up oc p 



